#####################################################################################################
# This script generates ideal points for North Carolina legislators. It estimates scores for each   #                                                                                                  #
# session of the state legislature. These scores are not in the same space (though the House and    #
# Senate for a given session are in the same space). Analyses across sessions requires              #
# inclusion of a session dummy.                                                                     #
#####################################################################################################

### Start by loading packages pscl package 
rm(list=ls(all=TRUE))
library(pscl)	

### Set directory
setwd("../Dataverse/Data/")

post1<-list()
post2<-list()
rc<-vector("list", 4)
### Set session year. To generate each session you need to change this (starts at 1846-47 and goes to 1858-59).
year<-"1858-59"

desc<- paste("NC ", year, "Both Chambers", sep=" ")
infile<-paste("matrices/Joint", year,".dat", sep="")
### Load roll call matrices using readKH (this assumes congress but works here as well)
rc <- readKH(infile)

### Drop roll calls that have fewer than 3 votes on one side, and drop legislators with fewer than 3 votes
rc.drop<-dropRollCall(rc,dropList=list(lop=3,legisMin=3))
summary(rc.drop)						#	Quick Summary
summary(rc.drop,verbose=TRUE)			#	Detailed Summary

### Estimate ideal points, normalized with a mean of 0 and SD of 1.
### 30,000 iterations, thinning every 30 iteration with 3000 burn in.
rc.ideal<-ideal(rc.drop,d=1,maxiter=30000,thin=30,burnin=3000,
                impute=FALSE,normalize=TRUE,store.item=TRUE,verbose=TRUE)
	

### Get the posterior means and distribution data and merge it with legislator info.
legdata.mean<-try(as.data.frame(summary(rc.ideal)$xm))
if(class(legdata.mean)=="try-error") { 
  xposterior<-rc.ideal$x	
  posteriorsd<-apply(xposterior,2,sd)
  legdata.mean<-rc.ideal$xbar
  legdata.sd<-as.data.frame(posteriorsd)
  posteriorHDR<-apply(xposterior,2,quantile, probs = c(.025, .975))
  post_vec<-t(posteriorHDR)
  legdata.HDR<-as.data.frame(post_vec)
} else {
  legdata.sd<-as.data.frame(summary(rc.ideal)$xsd)
  legdata.HDR<-as.data.frame(summary(rc.ideal)$xHDR)
}

new.ideal <- data.frame(rownames(legdata.mean), legdata.mean, legdata.sd, legdata.HDR)
leg.info<-rc.drop$legis
all.info<-cbind(leg.info, new.ideal )


plot(rc.ideal)

ideal.estimates<-all.info

### Save to a session specific file. These were then appended to each other to make "combined-session-by-session.csv"
outfile1<-paste("Joint", year, "IdealPoints.csv", sep="")
write.table(ideal.estimates, file = outfile1, sep = ",", col.names = NA)
